Pattern dynamics near homoclinic bifurcation in Rayleigh-Benard convection 
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We report for the first time the pattern dynamics in the vicinity of an inverse homoclinic bi- 
furcation in an extended dissipative system. We observe, in direct numerical simulations of three 
dimensional Rayleigh-Benard convection, a spontaneous breaking of a competition of two mutually 
perpendicular sets of oscillating cross rolls to one of two possible sets of oscillating cross rolls as 
the Rayleigh number is raised above a critical value. The time period of the cross-roll patterns 
diverges, and shows scaling behavior near the bifurcation point. This is an example of a transition 
from nonlocal to local pattern dynamics near an inverse homoclinic bifurcation. We also present a 
simple four-mode model that captures the pattern dynamics quite well. 



PACS numbers: 47.20.Ky, 47.55.pb, 47.20.Bp 



m 
o 

(N 
G 

00 
(N 



Oh 

d 



o 



x 



Extended dissipative systems driven away from ther- 
modynamic equilibrium often form patterns, if the driv- 
ing force exceeds a critical value Competing in- 
stabilities may lead to interesting pattern dynamics, 
which helps in understanding the underlying instabil- 
ity mechanism. Several patterns are observed in con- 
tinuum mechanical systems, such as Rayleigh-Benard 
systems @], Benard-Marangoni systems magneto- 
hydrodynamics J3] , ferrofluids @ , binary fluids @ , gran- 
ular materials 0| under shaking, biological systems Q , 
etc. Symmetries and dissipation play a very significant 
role in pattern selection in such systems [jj. The selec- 
tion of a pattern is a consequence of at least one broken 
symmetry of the system. Unbroken symmetries often in- 
troduce multiple patterns, which may lead to a transition 
from local to global pattern dynamics. The gluing 
of two limit cycles on two sides of a saddle point in the 
phase space of a given system is an example of a local to 
nonlocal bifurcation. It occurs when two limit cycles si- 
multaneously become homoclinic orbits of the same sad- 
dle point. This phenomenon has been recently observed 
in a variety of systems including liquid cryst als flll| , fluid 
dynamical systems [lU, biological systems jl3l|. optical 
systems (TJ, and electrical circuits [lU, and is a topic 
of current research. The pattern dynamics in the vicin- 
ity of a homoclinic bifurcation has, however, not been 
investigated in a fluid dynamical system. 

A Rayleigh-Benard system [l6|, [l?} , where a thin layer 
of a fluid is heated uniformly from below and cooled uni- 
formly from above, is a classical example of an extended 
dissipative system which shows a plethora of pattern- 
forming instabilities c haos fl8j . and turbulence fl9l ]. 
Low-Prandtl-number [20( and very low-Prandtl-numbcr 
convection [2l|, [22[ show three-dimensional oscillatory 
behavior close to the instability onset. In addition, 
the Rayleigh-Benard system possesses symmetries un- 
der translation and rotation in the horizontal plane that 
can introduce multiple sets of patterns. The possibil- 
ity of a homoclinic bifurcation and the pattern dynamics 
in its vicinity are unexplored in three dimensional (3D) 
Rayleigh-Benard convection. 



We report, in this article, for the first time the 
possibility of an inverse homoclinic bifurcation in di- 
rect numerical simulations (DNS) of three dimensional 
(3D) Rayleigh-Benard convection (RBC) in low-Prandtl- 
number fluids, and the results of our investigations of 
fluid patterns close to the bifurcation. We observe spon- 
taneous breaking of a periodic competition of two mutu- 
ally perpendicular sets of cross rolls to one set of oscillat- 
ing cross rolls, as the Rayleigh number Ra is raised above 
a critical value Rah- The time period of the oscillating 
patterns diverges, and shows scaling behavior in the close 
vicinity of the transition point. The exponents of scaling 
are asymmetric on the two sides of the transition point. 
We also present a simple four-mode model, which cap- 
tures not only the pattern dynamics in the vicinity of 
the inverse homoclinic bifurcation but also the whole se- 
quence of bifurcations observed in DNS quite well over 
a wide range of Rayleigh number in low-Prandtl-numbcr 
fluids. 

The hydrodynamics of RBC in a thin layer of Boussi- 
nesq fluid of thickness d, kinematic viscosity v, thermal 
diffusion coefficient k, and thermal expansion coefficient 
a, subjected to an adverse temperature gradient /?, is 
governed by the following dimensionless hydrodynamic 
equations: 



d t v + (v • V)v = 
Pr[d t e + 
V-v = 



-Vp + V z v + Ra6e 3 , 
(v-V)0] = V 2 6 
0, 



^3, 



(1) 
(2) 
(3) 



where v(x, y, z, t) = (ui,i>2> V3) is the velocity field, 
9(x,y,z,t) the convective temperature field, p the pres- 
sure due to convection, and e3 a unit vector directed 
against the direction of the acceleration due to gravity 
g. Lengths, time and temperature are measured in units 
of the fluid thickness d, viscous diffusion time d 2 /v, and 
v/3d/ ' k, respectively. We use thermally conducting and 
stress-free boundary conditions which imply that 6 = 
V3 = d z v\ = d z v 2 = at z = 0,1. All the fields are 
assumed to be periodic in the horizontal plane. The 
Rayleigh number Ra = a/3gd 4 jvK and Prandtl number 
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Pr = v / ti are two dimensionlcss numbers that decide 
the convective flow structures in the fluid. Convection 
appears when the reduced Rayleigh number r — Ra/Ra c 
with Ra c = 277r 4 /4 is raised above unity. 

We integrate the full hydrodynamic system (fT] - 
[3]) for low-Prandtl-numbcr (Pr < 0.025) fluids us- 
ing an object oriented code [23| based on pseudo- 
spectral method. The vertical velocity V3 and the 
temperature field 9 are expanded as: v 3 (x,y, z,t) = 
El, m ,n W l™n(t)e ik « x +™rt sin (nnz) and 6{x,y,z,t) = 

E;, m ,„ ®imn(t)e ik( - lx+my) sin (rnrz). The horizontal ve- 
locities v\ and V2 have similar expansion in xy plane with 
co-sinusoidal dependence in the z direction. The integers 
l,m,n can take values consistent with the equation of 
continuity (eq. [3]) and k = k c = tt / y/2. The size of the 
periodic cell for DNS is 2^2 x 2^2 x 1, and its resolu- 
tion is 64 x 64 x 64. The pattern dynamics is complex 
at the primary instability in very low-Prandtl- number 
fluids [22( due to chaotic flows just above the instabil- 
ity onset. We investigate pattern dynamics as soon as 
we observe the first oscillatory pattern close to the onset 
of convection. The simulation was started with random 
initial conditions, and it was continued until a steady 
state was reached. The steady state values of fields in 
a simulation were used as the initial conditions for the 
next simulation. The value of r was increased in small 
steps of size Ar (0.0001 < Ar < 0.01) and the numerical 
simulations were done for several values of r. We also 
repeated several runs starting with random initial condi- 
tions for different values of r and found no hysteresis in 
the parameter range considered here. Various observed 
convective patterns in the DNS for Pr = 0.01 and Pr = 
are listed in the first two columns of Table |U 



TABLE I: Convective patterns computed from DNS in two 
columns in the middle and from a model in the last column. 



Convective 
patterns 


DNS 


Model 


r(Pr = 0.01) 


r(Pr = 0) 


r(Pr = 0) 


OCR-I 


1.010 - 1.0835 


1.0049 - 1.0708 


1.010 - 1.0953 


OCR-II 


1.0836 - 1.1200 


1.0709 - 1.1315 


1.0954 - 1.1584 


CR 


1.1210 - 1.1990 


1.1316 - 1.2005 


1.1585 - 1.2519 


SQ 


1.2000 - 1.4200 


1.2006 - 1.4297 


1.2520 - 1.5128 



The first ordered state for Pr = 0.01 appears in the 
form of a competition of two mutually perpendicular sets 
of oscillatory cross rolls (OCR-I) which continues to exist 
until r — 1.0835. The first row of Fig. [TJ displays the pat- 
tern dynamics for r = 1.076. Two sets of cross rolls, one 
with |Wioi| > I Won I and another with |Wi i| < |Won|, 
appear periodically. The patterns appear as squares 
when the amplitudes of the two sets of cross rolls be- 
come equal. The competition represents a global pattern 
dynamics. As r is raised above a critical r^ = 1.0835, the 
competing cross rolls (OCR-I) spontaneously break into 
two possible oscillatory cross rolls (OCR-II). The second 




FIG. 1: Contour plots of the temperature field at the mid 
plane (z = 0.5) near inverse homoclinic bifurcation as ob- 
served in DNS (Pr = 0.01). The upper row (r = 1.076) 
shows competition between two sets of oscillating cross rolls. 
The middle and lower rows show two possibilities of oscillating 
cross rolls (r — 1.088). 



and third rows of Fig. [TJ show two possibilities of OCR-II 
for r = 1.088. We get oscillating cross-roll patterns with 
either |Wi i| > | Won I or W101 < Won. Two sets of mul- 
tiple solutions, which are connected by rotation about a 
vertical axis by 7r/2, continue until r = 1.120. Further 
increase in r leads to the appearance of two sets of sta- 
tionary cross rolls (CR) at r = 1.121, which is observed 
till r = 1.199. Raising the value of r even further leads 
to a transition from stationary cross rolls (CR) to sta- 
tionary square (SQ) patterns. The similar sequence is 
observed in the limit of Pr — > (sec Table [TJ) . The range 
of competing cross rolls becomes wider as Pr decreases. 
We have observed the spontaneous breaking of compet- 
ing cross rolls to two sets of oscillatory cross rolls in fluids 
with < Pr < 0.025. 

Figure [2] shows the details of a transition from global 
to local pattern dynamics for Pr = 0.01. The divergence 
of the time period of oscillatory patterns close to the 
transition point is displayed in Fig. [2Ji. The amplitude 
of the largest Fourier mode of OCR-I patterns decreases 
linearly with the increase of r and shows two possible 
values just above the transition (r = r/j) point (Fig. [2b). 
The appearance of two amplitudes signifies a transition 
from a nonlocal to a local pattern dynamics. The scaling 
of time period r of oscillating pattern on both sides of the 
transition point is displayed in Fig.[5J:, showing asymme- 
try. The time period r of competing patterns scales with 
e = \r — r/j| as e~° 115 before transition and as £ -° 033 
after transition. The scaling behavior of the time period 
of OCR sug gest s the transition to be inverse homoclinic. 
Simulations [2l| of low- Pr RBC with no-slip boundaries 
are known to show relaxation oscillation of patterns. The 
homoclinic instability may be accessible to experiments, 
if performed in square containers, in this regime. 

We now construct a simple low dimensional model to 
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FIG. 2: Scaling behavior near inverse homoclinic bifurcation 
(Pr = 0.01) as obtained from DNS: (a) The divergence of the 
dimesionless time period r of OCR patterns close to the bi- 
furcation, (b) The variation of maxima of Fourier mode Wioi 
with the reduced Rayleigh number r near the bifurcation, 
showing the spontaneous transition from a global oscillation 
(star *) to two possible local oscillations(cross x) at the ho- 
moclinic point, (c) The time period r of OCR patterns scales 
with e = \r — | . The scaling exponents are different before 
and after the bifurcation. 



analyze pattern dynamics near the inverse homoclinic bi- 
furcation. For this purpose we take the limit of vanishing 
Prandtl number (Pr — > 0). As the temperature field is 
slaved to the vertical velocity, the number of modes rep- 
resenting the effective dynamics is expected to be smaller 
in this limit. We begin with the standard Galerkin tech- 
nique to derive a low-dimensional-model [22j. We ex- 
pand the vertical velocity v 3 and the vertical vorticity 
Z = (V x v)-e 3 such that the essential modes to describe 
two sets of mutually perpendicular rolls, cross rolls, and 
the nonlinear interaction between them are retained. We 
keep five velocity modes Wioi, Won, W211, W121, and 
W112, and two vorticity modes Zhq and Zn 2 . The hy- 
drodynamic equations are projected on these modes. We 
then adiabatically eliminate modes W112, -Zaio and Zn 2 . 
This leads to a simple four-mode model given by, 

X = fi 1 X + X 1 X 2 A(a 1 X + a 2 Y)+Y 1 Y 2 A(a 3 X + a 4 Y) 
+ a 5 [XlY ll X 2 l Y 2 ] T + a^X x Y^X 2 Yl\ J ', 

Y = /i 2 Y + X 1 X 2 A(b 1 X + b 2 Y)+Y 1 Y 2 A(b 3 X + b 4 Y) 
+ b 5 [XiY 1 ,XfY 2 ] T + b 6 [X 1 Yi,X 2 Y 1 2 ] T , (4) 

where X = [X U X 2 ] T = [W 1Q i, W 011 ] T , Y = [Y X ,Y 2 ] T = 
[W 121 ,W 211 ] T , A = [0 1;1 0], /ii = 3^ 2 (r - l)/2 
and fi 2 = 7r 2 (135r — 343)/98. The coefficients are: 
ai = -3/100, a 2 = 31/3000, a 3 = -209/30000, a 4 = 
63/60000, a 5 = -47/1500, a 6 = 1/200, h = -93/700, 
b 2 = -67/7000, 63 = -7407/70000, b A = -3969/700000, 
65 = -928/7000, and b 6 = 3816/70000. The superscript 
T denotes transpose of a matrix. The model is valid for 
r < 343/135 (i.e., fj, 2 < 0). 

The model is integrated using standard fourth order 
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FIG. 3: Comparison of results of the model with those from 
DNS for Pr — 0: Scaling of the dimensionless time period 
r of the oscillating patterns with e = \r — ru\ before (solid 
curve) and after (dashed curve) the bifurcation, as observed 
in DNS (a) and the model (b). The exponents from DNS 
are: 71 = 0.086 and 72 = 0.158, while those from the model 
are: 71 = 0.092 and 72 = 0.124. Temporal variation of the 
Fourier mode W101 computed from DNS before (c) and after 
(d) the homoclinic bifurcation. Temporal variation of the 
mode Wwi computed from the model before (e) and after (f) 
the bifurcation. 



Runge Kutta (RK4) method. The second and third 
columns of Table U summarize the results obtained from 
DNS for Pr = and the model, respectively. The model 
captures the sequence of bifurcations quite accurately in 
a wide range of r as observed in DNS. The difference 
in the lower and the upper bounds for the range of r 
for any solution computed from the model and DNS is 
within 6%. 

Figure [3] gives the comparison of results obtained from 
the model and those from DNS for Pr = near homo- 
clinic bifurcation. The time period r of the competing 
cross rolls scales with e = |r — as e~ 71 before the 
transition (solid curve) and as e~ 72 after the transition 
(dashed curve). The values of the scaling exponent 71 ob- 
tained from the DNS (Fig. |3ja,) and the model (Fig. [3b) 
are 0.086 and 0.092 respectively. Two exponents are dif- 
ferent showing asymmetry in scaling behavior. DNS and 
the model yield 72 equal to 0.158 and 0.124 respectively. 
Figurcs[3j; & d show the temporal variation of the Fourier 
mode W101 obtained from DNS before and after the bi- 
furcation. The similar behavior is observed in the model 
(fig- 12k, f)- The model reveals that the unstable SQ pat- 
terns exist as saddle fixed points for 1 < r < 1.252 but 
become stable at r = 1.252. The competing cross rolls 
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(OCR-I) break into two possible sets of OCR-II when the 
amplitude of OCR-I oscillation touches a saddle square. 
This confirms the transition from OCR-I to OCR-II as 
an inverse homoclinic bifurcation. The model includes 
very few modes, and therefore shows higher values of 
the Fourier mode Wioi- The time periods of oscillating 
patterns obtained from the model and DNS are in good 
agreement both before and after the transition. 

We have investigated the possibility of an inverse ho- 
moclinic bifurcation and associated pattern dynamics in 
RBC. The spontaneous breaking of a competition be- 
tween two mutually perpendicular sets of oscillatory cross 



rolls into one set of oscillatory cross rolls occurs close 
to the bifurcation point. The time period of oscillation 
shows asymmetric scaling on the two sides of the bifurca- 
tion point. We have also constructed a simple four-mode 
model which captures accurately the sequence of bifurca- 
tions including the pattern dynamics near the homoclinic 
bifurcation. The model with different values of the coeffi- 
cients cii and bi may be useful to study pattern dynamics 
on square lattices in other extended systems with similar 
symmetries. 

We have benefitted from fruitful discussions with J.K. 
Bhattacharjee, H. Pharasi, L.K. Dey, and D. Kumar. 
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